################################################################################
################################################################################
## Uganda National Panel Survey (UNPS)
## UNPS and Aid data maps
## Ray Caraher
## Date started: 2021-05-18
################################################################################
################################################################################

print(paste(Sys.time(), Sys.getenv("USERNAME")))

library(sf)
library(sp)
library(geosphere)
library(rnaturalearth)
library(rnaturalearthdata)
library(rgeos)
library(mapview)
library(rgdal)
library(ggrepel)
library(ggmap)
library(broom)
library(tidyverse)
library(leaflet)
library(leaflet.extras)
library(htmltools)
library(ggmap)
library(cowplot)
library(openxlsx)
library(rio)
library(openxlsx)
library(tidyverse)
library(tidylog)
library(stargazer)

options(scipen = 999)

rm(list = ls())

sys_info <- paste(Sys.info()[["sysname"]], 
                  Sys.info()[["version"]], 
                  Sys.info()[["user"]], sep = "-")

################################################################################
## Read in merged UNPS data

## Set base directory to relative location to read file

base_directory <- "C:/Users/raymo/OneDrive - University of Massachusetts/Academic/RA/UNPS/Analysis/"

data_directory <- paste0(base_directory, "replication-data/")

## Import GSEC data

filename <- paste0(data_directory, "unps_gsec_aid_rep_data.dta")

gsec_aid <- import(filename)

## Import data on aid by district

filename <- paste0(data_directory, "aid_by_districts.csv")

aid_dist <- import(filename)

aid_dist <- aid_dist %>%
  mutate(district = ifelse(district == "SEMBABULE", "SSEMBABULE", district))

  ## Import aid projects by coordinates

filename <- paste0(data_directory, "aid_by_coords.csv")
aid_coords <- import(filename)

## Filter to WASH aid

ws_coords <- aid_coords %>%
  filter(wat_san == 1)

ws_coords_pc1 <- filter(ws_coords, precision_code == 1)

################################################################################
## Average, min, max distance to projects across years

# dist_cols <- c("min_dist_2009_proj", "min_dist_2010_proj", "min_dist_2011_proj",
#                "min_dist_2012_proj", "min_dist_2013_proj", "min_dist_2014_proj")

# gsec_aid <- gsec_aid %>%
#   rowwise() %>%
#   mutate(min_min_dist = min(c_across(all_of(dist_cols))),
#          max_min_dist = max(c_across(all_of(dist_cols))),
#          mean_min_dist = mean(c_across(all_of(dist_cols)))
#            )

# dists_by_region <- gsec_aid %>%
#   group_by(region) %>%
#   select(min_min_dist, max_min_dist, mean_min_dist) %>%
#   summarise_all(mean, na.rm = T)


################################################################################

## Focus on 2015 data

gsec_aid_1516 <- filter(gsec_aid, year == "2015-2016")

## Read in shape file

directory <- paste0(data_directory, "uga_admbnda_ubos_20200824_shp/")

filename <- paste0(directory, "uga_admbnda_adm2_ubos_20200824.shp")

uga_dists <- st_read(filename)

gsec_aid_1516 <- gsec_aid_1516 %>%
  mutate(label = id,
         district = toupper(district))

uga_dists_df <- readOGR(filename)
uga_dists_df <- tidy(uga_dists_df, region = "ADM2_EN")
uga_dists_df <- uga_dists_df %>%
  mutate(district = toupper(id))


################################################################################
## Distribution of HH across districts

## Plot districts

p1 <- ggplot(uga_dists_df) +
  geom_polygon(aes(x = long, y = lat, group = id), 
               color = "black", fill = "lightgray") +
  theme_void() +
  coord_map()

print(p1)

gsec_hh <- gsec_aid_1516 %>%
  distinct(id, .keep_all = T)

gsec_hh_dist <- gsec_hh %>%
  group_by(district) %>%
  summarise(num_hh = n())

match(gsec_hh_dist$district, uga_dists_df$district)

uga_dists_df <- left_join(uga_dists_df,
                          gsec_hh_dist, by = "district")

  
p1 <- ggplot(uga_dists_df) +
  geom_polygon(aes(x = long, y = lat, group = id, fill = num_hh), 
               color = "black") +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Households") +
  theme_void() +
  coord_map()

print(p1)

## Aid by district

uga_dists_df <- left_join(uga_dists_df,
                          aid_dist, by = "district")

p2 <- ggplot(uga_dists_df) +
  geom_polygon(aes(x = long, y = lat, group = id, fill = `code_1_Water and sanitation`), 
               color = "black") +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Aid projects") +
  theme_void() +
  coord_map()

print(p2)

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

ws_coords_pc1.sf <- st_as_sf(ws_coords_pc1, coords = c("longitude", "latitude"), 
                          crs = 4326, agr = "constant")


uga_dists <- uga_dists %>%
  mutate(district = toupper(ADM2_EN))

uga_dists <- left_join(uga_dists,
                          gsec_hh_dist, by = "district")

uga_dists.sf <- st_as_sf(uga_dists)

## Plot aid and households

p3 <- ggplot(data = world) +
  geom_sf() +
  geom_sf(data = uga_dists.sf) +
  geom_sf(data = ws_coords_pc1.sf) +
  coord_sf(xlim = c(28, 36), ylim = c(-1.5, 4.5), expand = F) +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Aid projects")

print(p3)

p4 <- ggplot(data = world) +
  geom_sf() +
  geom_sf(data = uga_dists.sf, aes(fill = num_hh)) +
  geom_sf(data = ws_coords_pc1.sf, size = 2, shape = 17) +
  coord_sf(xlim = c(28, 36), ylim = c(-1.5, 4.5), expand = F) +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Households")

print(p4)

## Add labels

text <- data.frame(place = c("Kampala", "D.R. Congo", "South Sudan", "Kenya", "Tanzania"),
                   lat = c(0.5, 2, 4, 0, -1.3),
                   long = c(33, 30, 32, 34.5, 32.5))

text.sf <- st_as_sf(text, coords = c("long", "lat"), 
                    crs = 4326, agr = "constant", remove = F)

p5 <- ggplot(data = world) +
  geom_sf() +
  geom_sf(data = uga_dists.sf, aes(fill = num_hh)) +
  geom_sf(data = ws_coords_pc1.sf, size = 2, shape = 17) +
  geom_text(data = text.sf, aes(x = long, y = lat, label = place), 
            size = 2.5, col = "black", fontface = "bold") +
  coord_sf(xlim = c(29, 35.5), ylim = c(-1.5, 4.5), expand = F) +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Households") +
  labs(x = "Longitude", y = "Latitude")

print(p5)

addSmallLegend <- function(myPlot, pointSize = 4, textSize = 8, spaceLegend = 1) {
  myPlot +
    guides(shape = guide_legend(override.aes = list(size = pointSize)),
           color = guide_legend(override.aes = list(size = pointSize))) +
    theme(legend.title = element_text(size = textSize + 2), 
          legend.text  = element_text(size = textSize),
          legend.key.size = unit(spaceLegend, "lines"))
}

p6 <- addSmallLegend(p5)

print(p6)

## All aid projects

aid_coords_pc1 <- filter(aid_coords, precision_code == 1)
aid_coords_pc1 <- aid_coords_pc1 %>%
  mutate(type2 = ifelse(type == "Agriculture", "Agriculture, Logisitics, and Other",
                 ifelse(type == "Logistics", "Agriculture, Logisitics, and Other",
                 ifelse(type == "Other", "Agriculture, Logisitics, and Other", type))),
         type2.f = as_factor(type2))

aid_coords_pc1.sf <- st_as_sf(aid_coords_pc1, coords = c("longitude", "latitude"), 
                             crs = 4326, agr = "constant")



p7 <- ggplot(data = world) +
  geom_sf() +
  geom_sf(data = uga_dists.sf, aes(fill = num_hh)) +
  geom_sf(data = aid_coords_pc1.sf, size = 2, aes(color = type2.f)) +
  geom_text(data = text.sf, aes(x = long, y = lat, label = place), 
            size = 2.5, col = "black", fontface = "bold") +
  coord_sf(xlim = c(29, 35.5), ylim = c(-1.5, 4.5), expand = F) +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Households") +
  scale_color_discrete(name = "Aid type") +
  labs(x = "Longitude", y = "Latitude")

p7 <- addSmallLegend(p7)

p7 <- p7 + theme(legend.position = "bottom")

print(p7)

## Just Water and Sanitation Aid

aid_coords_pc1_ws <- filter(aid_coords, precision_code == 1, type == "Water and sanitation")
aid_coords_pc1_ws.sf <- st_as_sf(aid_coords_pc1_ws, coords = c("longitude", "latitude"), 
                                 crs = 4326, agr = "constant")

aid_coords_pc1_proj <- aid_coords_pc1_ws %>% group_by(project_title) %>%
  summarise(proj_count = n(),
            proj_year_end = mean(transactions_end_year),
            comm = mean(total_commitments),
            disp = mean(total_disbursements),
            code = mean(precision_code)
  ) %>%
  ungroup()

p7b <- ggplot(data = world) +
  geom_sf() +
  geom_sf(data = uga_dists.sf, aes(fill = num_hh)) +
  geom_sf(data = aid_coords_pc1_ws.sf, size = 2, color = "blue") +
  geom_text(data = text.sf, aes(x = long, y = lat, label = place), 
            size = 2.5, col = "black", fontface = "bold") +
  coord_sf(xlim = c(29, 35.5), ylim = c(-1.5, 4.5), expand = F) +
  scale_fill_gradient(low = "white", high = "black",
                      name = "Households") +
  scale_color_discrete(name = "Aid type") +
  labs(x = "Longitude", y = "Latitude")

p7b <- addSmallLegend(p7b)

p7b <- p7b + theme(legend.position = "bottom")

print(p7b)


p8 <- ggplot(data = aid_coords_pc1_ws) +
  geom_bar(aes(x = transactions_end_year)) +
scale_x_continuous(breaks = seq(2009, 2014)) +
  scale_y_continuous(breaks = seq(0, 20, 2)) +
  labs(x = "Aid project completion", y = "Count")
  
print(p8)

## Figure of treatment status

## Transform to long

treatments <- gsec_aid %>%
  select(id, pid, year, treat_1km, treat_5km, treat_10km, treat_15km, treat_20km, treat_25km,
         treatment_year_1km, treatment_year_5km, treatment_year_10km, treatment_year_15km,
         treatment_year_20km, treatment_year_25km)

#pdf(paste0(base_directory, "/../figures/treatment-status-bargraphs.pdf"))

pp1 <- ggplot() +
  geom_bar(aes(x = treat_25km, y = year, 
               fill = as.factor(treatment_year_25km)),
           stat = "identity",
           data = filter(gsec_aid, treat_25km == 1)) +
  labs(y = "Survey wave", x = "Treated - 25km", fill = "Treatment year") +
  theme(legend.position = "bottom",
        text = element_text(size = 7))

legend <- get_legend(pp1 +
                     guides(color = guide_legend(nrow = 1)))

pp1 <- pp1 + theme(legend.position = "none")

print(pp1)

pp2 <- ggplot() +
  geom_bar(aes(x = treat_20km, y = year, 
               group = treatment_year_20km, fill = as.factor(treatment_year_20km)),
           stat = "identity",
           data = filter(gsec_aid, treat_20km == 1)) +
  labs(y = "Survey wave", x = "Treated - 20km", fill = "Treatment year")+
  theme(legend.position = "none",
        text = element_text(size = 7))

print(pp2)

pp3 <- ggplot() +
  geom_bar(aes(x = treat_15km, y = year, 
               group = treatment_year_15km, fill = as.factor(treatment_year_15km)),
           stat = "identity",
           data = filter(gsec_aid, treat_15km == 1)) +
  labs(y = "Survey wave", x = "Treated - 15km", fill = "Treatment year")+
  theme(legend.position = "none",
        text = element_text(size = 7))

print(pp3)

pp4 <- ggplot() +
  geom_bar(aes(x = treat_10km, y = year, 
               group = treatment_year_10km, fill = as.factor(treatment_year_10km)),
           stat = "identity",
           data = filter(gsec_aid, treat_10km == 1)) +
  labs(y = "Survey wave", x = "Treated - 10km", fill = "Treatment year")+
  theme(legend.position = "none",
        text = element_text(size = 7))

print(pp4)

pp5 <- ggplot() +
  geom_bar(aes(x = treat_5km, y = year, 
               group = treatment_year_5km, fill = as.factor(treatment_year_5km)),
           stat = "identity",
           data = filter(gsec_aid, treat_5km == 1)) +
  labs(y = "Survey wave", x = "Treated - 5km", fill = "Treatment year")+
  theme(legend.position = "none",
        text = element_text(size = 7))

print(pp5)

pp6 <- ggplot() +
  geom_bar(aes(x = treat_1km, y = year, 
               group = treatment_year_1km, fill = as.factor(treatment_year_1km)),
           stat = "identity",
           data = filter(gsec_aid, treat_1km == 1)) +
  labs(y = "Survey wave", x = "Treated - 1km", fill = "Treatment year")+
  theme(legend.position = "none",
        text = element_text(size = 7))

print(pp6)

pg1 <- plot_grid(pp6,pp5,pp4,pp3,pp2,pp1, nrow = 3, ncol = 2)
print(pg1)
pg <- plot_grid(pg1, legend, nrow = 2, 
                rel_widths = c(1, 0.25), rel_heights = c(1, 0.25))
print(pg)


################################################################################
## Export Figures

figures_directory <- paste0(base_directory, "replication-figures/")

save_plot(p7b, 
          base_height = 6,
          filename = paste0(figures_directory, "map_ws_hh.png"))

save_plot(p8,
          filename = paste0(figures_directory, "aid_project_completion.png"))

save_plot(pg,
          filename = paste0(figures_directory, "treatment-status-bargraphs.png"))


################################################################################
## Misc. tables
################################################################################
## Set tables directory

tables_directory <- paste0(base_directory, "replication-tables/")

## Table 1 - Average aid project disbursements

avg_disb <- aid_coords_pc1_ws %>%
  group_by(transactions_end_year) %>%
  summarise(count_proj = n(),
    mean_tot_disb = round(mean(total_disbursements),1),
    meant_tot_comm = round(mean(total_commitments),1)
)


colnames(avg_disb) <- c("Project end year", "Number of projects", "Avg. disbursements", "Avg. commitments")


stargazer(avg_disb, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Projects by year and funding",
          label = "proj_count_funds",
          out = paste0(tables_directory, "ws_project_count_funds.tex"))

stargazer(avg_disb, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Projects by year and funding",
          label = "proj_count_funds",
          type = "html",
          out = paste0(tables_directory, "ws_project_count_funds.rtf"))

################################################################################
## Create variable table

filename <- paste0(data_directory, "variables.xlsx")

var_tab <- import(filename)

stargazer(var_tab, summary = F,
          font.size = "small",
          column.sep.width = "1pt",
          title = "Variables",
          label = "var_tab",
          out = paste0(tables_directory, "var_tab.tex"))

################################################################################
## Create table for average distances - not available due to data confidentiality reasons

# dists_by_region <- dists_by_region %>%
#   filter(!is.na(region)) %>%
#   mutate(region = ifelse(region == 0, "Kampala",
#                   ifelse(region == 1, "Central",
#                   ifelse(region == 2, "Eastern",
#                   ifelse(region == 3, "Northern",
#                   ifelse(region == 4, "Western", "ERR")))))) %>%
#   mutate_if(is.numeric, ~round(., 1))

# colnames(dists_by_region) <- c("Region", "Avg. min. distance to closest project",
#                                "Max. distance to closest project",
#                                "Avg. distance to closest project")

# dists_by_region <- select(dists_by_region, Region, `Avg. min. distance to closest project`)

# stargazer(dists_by_region, summary = F,
#           font.size = "small",
#           column.sep.width = "2pt",
#           title = "Avg. min. distance to closest project by region",
#           label = "dist_to_proj_reg",
#           out = paste0(tables_directory, "dist_to_proj_reg.tex"))

# stargazer(dists_by_region, summary = F,
#           font.size = "small",
#           column.sep.width = "2pt",
#           title = "Avg. min. distance to closest project by region",
#           label = "dist_to_proj_reg",
#           type = "html",
#           out = paste0(tables_directory, "dist_to_proj_reg.rtf"))

################################################################################
## Create table for access to water

## Improved water

water_improved_tab <- gsec_aid %>%
  group_by(region, year) %>%
  summarise(water_imp_mean = mean(water_improved, na.rm = T))

water_improved_tab <- water_improved_tab %>% filter(!is.na(region)) %>%
  mutate(region = ifelse(region == 0, "Kampala",
                  ifelse(region == 1, "Central",
                  ifelse(region == 2, "Eastern",
                  ifelse(region == 3, "Northern",
                  ifelse(region == 4, "Western", "ERR")))))) %>%
  mutate_if(is.numeric, ~round(., 3))

water_improved_tab <- water_improved_tab %>%
  pivot_wider(names_from = region, values_from = water_imp_mean)

gsec_aid <- gsec_aid %>%
  group_by(year, id) %>%
  mutate(hh_year_id = cur_group_id()) %>%
  ungroup()

gsec_aid_hh <- gsec_aid %>%
  distinct(hh_year_id, .keep_all = T)

water_improved_tab_hh <- gsec_aid_hh %>%
  group_by(region, year) %>%
  summarise(water_imp_mean = mean(water_improved, na.rm = T))


water_improved_tab_hh <- water_improved_tab_hh %>% filter(!is.na(region)) %>%
  mutate(region = ifelse(region == 0, "Kampala",
                         ifelse(region == 1, "Central",
                                ifelse(region == 2, "Eastern",
                                       ifelse(region == 3, "Northern",
                                              ifelse(region == 4, "Western", "ERR")))))) %>%
  mutate_if(is.numeric, ~round(., 3))

water_improved_tab_hh <- water_improved_tab_hh %>%
  pivot_wider(names_from = region, values_from = water_imp_mean)

colnames(water_improved_tab_hh) <- c("Year", colnames(water_improved_tab_hh)[-1])

water_improved_tab_hh <- water_improved_tab_hh %>%
  mutate_if(is.numeric, ~round(., 1))

## Water waiting time

water_wait_tab <- gsec_aid %>%
  group_by(region, year) %>%
  summarise(water_wait_mean = mean(water_waiting_time, na.rm = T))

mean(gsec_aid$water_waiting_time, na.rm = T)
mean(gsec_aid_hh$water_waiting_time, na.rm = T)

water_wait_tab <- water_wait_tab %>% filter(!is.na(region)) %>%
  mutate(region = ifelse(region == 0, "Kampala",
                         ifelse(region == 1, "Central",
                                ifelse(region == 2, "Eastern",
                                       ifelse(region == 3, "Northern",
                                              ifelse(region == 4, "Western", "ERR")))))) %>%
  mutate_if(is.numeric, ~round(., 3))

water_wait_tab <- water_wait_tab %>%
  pivot_wider(names_from = region, values_from = water_wait_mean)

water_wait_tab_hh <- gsec_aid_hh %>%
  group_by(region, year) %>%
  summarise(water_wait_mean = mean(water_waiting_time, na.rm = T))


water_wait_tab_hh <- water_wait_tab_hh %>% filter(!is.na(region)) %>%
  mutate(region = ifelse(region == 0, "Kampala",
                         ifelse(region == 1, "Central",
                                ifelse(region == 2, "Eastern",
                                       ifelse(region == 3, "Northern",
                                              ifelse(region == 4, "Western", "ERR")))))) %>%
  mutate_if(is.numeric, ~round(., 3))

water_wait_tab_hh <- water_wait_tab_hh %>%
  pivot_wider(names_from = region, values_from = water_wait_mean)

colnames(water_wait_tab_hh) <- c("Year", colnames(water_wait_tab_hh)[-1])

water_wait_tab_hh <- water_wait_tab_hh %>%
  mutate_if(is.numeric, ~round(., 1))


## Water travel time

water_trav_tab <- gsec_aid %>%
  group_by(region, year) %>%
  summarise(water_trav_mean = mean(water_travel_time, na.rm = T))

mean(gsec_aid$water_travel_time, na.rm = T)
mean(gsec_aid_hh$water_travel_time, na.rm = T)

water_trav_tab <- water_trav_tab %>% filter(!is.na(region)) %>%
  mutate(region = ifelse(region == 0, "Kampala",
                         ifelse(region == 1, "Central",
                                ifelse(region == 2, "Eastern",
                                       ifelse(region == 3, "Northern",
                                              ifelse(region == 4, "Western", "ERR")))))) %>%
  mutate_if(is.numeric, ~round(., 3))

water_trav_tab <- water_trav_tab %>%
  pivot_wider(names_from = region, values_from = water_trav_mean)

water_trav_tab_hh <- gsec_aid_hh %>%
  group_by(region, year) %>%
  summarise(water_trav_mean = mean(water_travel_time, na.rm = T))


water_trav_tab_hh <- water_trav_tab_hh %>% filter(!is.na(region)) %>%
  mutate(region = ifelse(region == 0, "Kampala",
                         ifelse(region == 1, "Central",
                                ifelse(region == 2, "Eastern",
                                       ifelse(region == 3, "Northern",
                                              ifelse(region == 4, "Western", "ERR")))))) %>%
  mutate_if(is.numeric, ~round(., 3))

water_trav_tab_hh <- water_trav_tab_hh %>%
  pivot_wider(names_from = region, values_from = water_trav_mean)

colnames(water_trav_tab_hh) <- c("Year", colnames(water_wait_tab_hh)[-1])

water_trav_tab_hh <- water_trav_tab_hh %>%
  mutate_if(is.numeric, ~round(., 1))

################################################################################
## Export tables

stargazer(water_improved_tab_hh, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Proportion of households with improved water, by region and survey wave",
          label = "water_improved_tab_reg_year",
          out = paste0(tables_directory, "water_improved_tab_reg_year.tex"))

stargazer(water_improved_tab_hh, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Proportion of households with improved water, by region and survey wave",
          label = "water_improved_tab_reg_year",
          type = "html",
          out = paste0(tables_directory, "water_improved_tab_reg_year.rtf"))

stargazer(water_trav_tab_hh, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Mean household travel time to water, by region and survey wave",
          label = "water_trav_tab_reg_year",
          out = paste0(tables_directory, "water_trav_tab_reg_year.tex"))

stargazer(water_trav_tab_hh, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Mean household travel time to water, by region and survey wave",
          label = "water_trav_tab_reg_year",
          type = "html",
          out = paste0(tables_directory, "water_trav_tab_reg_year.rtf"))

stargazer(water_wait_tab_hh, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Mean household waiting time for water, by region and survey wave",
          label = "water_wait_tab_reg_year",
          out = paste0(tables_directory, "water_wait_tab_reg_year.tex"))

stargazer(water_wait_tab_hh, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          title = "Mean household waiting time for water, by region and survey wave",
          label = "water_wait_tab_reg_year",
          type = "html",
          out = paste0(tables_directory, "water_wait_tab_reg_year.rtf"))

################################################################################
## Examine trends in ODA data for Uganda

## Import data

filename <- paste0(data_directory, "crs_uga_sect_90_20.csv")

crs_uga_sect_90_20 <- import(filename)

## Find total aid

crs_uga_tot_90_20 <- crs_uga_sect_90_20 %>%
    group_by(Year, RecipientCode, RecipientName) %>%
    summarise(All_USD_Commitment = sum(USD_Commitment, na.rm = T),
        All_USD_Disbursement = sum(USD_Disbursement, na.rm = T),
        All_USD_Commitment_Defl = sum(USD_Commitment_Defl, na.rm = T),
        All_USD_Disbursement_Defl = sum(USD_Disbursement_Defl, na.rm = T))

## Find WASH aid

crs_uga_wash_90_20 <- crs_uga_sect_90_20 %>%
    filter(SectorCode == 140) %>%
    rename(WASH_USD_Commitment = USD_Commitment,
        WASH_USD_Disbursement = USD_Disbursement,
        WASH_USD_Commitment_Defl = USD_Commitment_Defl,
        WASH_USD_Disbursement_Defl = USD_Disbursement_Defl
        ) %>%
    select(-c("SectorCode", "SectorName"))

## Combine

crs_uga_oda <- crs_uga_tot_90_20 %>%
  left_join(crs_uga_wash_90_20,
    by = c("Year", "RecipientCode", "RecipientName"))

## Create variables

crs_uga_oda <- crs_uga_oda %>%
  ungroup()

crs_uga_oda <- crs_uga_oda %>%
  mutate(wash_commit_prop = WASH_USD_Commitment_Defl/All_USD_Commitment_Defl,
    wash_disb_prop = WASH_USD_Disbursement_Defl/All_USD_Disbursement_Defl,
    wash_commit_pc = WASH_USD_Commitment_Defl/pop,
    wash_disb_pc = WASH_USD_Disbursement_Defl/pop)

## Read in Basic Water Supply data

filename <- paste0(data_directory, "crs_uga_bw_10_20.csv")

crs_uga_bw_10_20 <- import(filename)

crs_uga_oda <- crs_uga_oda %>%
  left_join(crs_uga_bw_10_20, by = "Year")

crs_uga_oda <- crs_uga_oda %>%
  mutate(bw_disb_prop = `BW Disbursments`/WASH_USD_Disbursement_Defl,
    bw_commit_prop = `BW Commitments`/WASH_USD_Commitment_Defl,
    bw_disb_pc = `BW Disbursments`/pop,
    bw_commit_pc = `BW Commitments`/pop)

## Collapse to 5-year periods

crs_uga_oda <- crs_uga_oda %>%
  arrange(Year) %>%
  mutate(
    period = gl(ceiling(n() / 5), 5, length = n()),
    period = case_when(period == 1 ~ "1990-94",
      period == 2 ~ "1995-99",
      period == 3 ~ "2000-04",
      period == 4 ~ "2005-09",
      period == 5 ~ "2010-14",
      period == 6 ~ "2015-20",
      period == 7 ~ "2015-20"))

## Make table

crs_uga_oda_5yr <- crs_uga_oda %>%
  select(-c("Year", "RecipientCode", "RecipientName")) %>%
  group_by(period) %>%
  summarise_all(mean, na.rm = T)

crs_uga_oda_5yr <- crs_uga_oda_5yr %>%
  rename("Avg. Total ODA Commitments" = "All_USD_Commitment_Defl",
  "Avg. Total ODA Disbursements" = "All_USD_Disbursement_Defl",
  "Avg. WASH ODA Commitments" = "WASH_USD_Commitment_Defl",
  "Avg. WASH ODA Disbursements" = "WASH_USD_Disbursement_Defl",
  "Avg. WASH Commitments Share" = "wash_commit_prop",
  "Avg. WASH Disbursements Share" = "wash_disb_prop",
  "Avg. WASH Commitments per capita" = "wash_commit_pc",
  "Avg. WASH Disbursements per capita" = "wash_disb_pc",
  "Avg. BW Disbursements" = "BW Disbursments",
  "Avg. BW Commitments" = "BW Commitments",
  "Avg. BW Disbursements Share of WASH" = "bw_disb_prop",
  "Avg. BW Commitments Share of WASH" = "bw_commit_prop",
  "Avg. BW Disbursements per capita" = "bw_disb_pc",
  "Avg. BW Commitments per capita" = "bw_commit_pc",
  "Period" = "period"
  ) %>%
  select(Period,
  "Avg. Total ODA Commitments", "Avg. Total ODA Disbursements",
  "Avg. WASH ODA Commitments","Avg. WASH ODA Disbursements", 
  "Avg. WASH Commitments Share", "Avg. WASH Disbursements Share",
  "Avg. WASH Commitments per capita", "Avg. WASH Disbursements per capita",
  "Avg. BW Commitments", "Avg. BW Disbursements",
  "Avg. BW Commitments Share of WASH", "Avg. BW Disbursements Share of WASH",
  "Avg. BW Commitments per capita", "Avg. BW Disbursements per capita"
  )

crs_uga_oda_5yr <- crs_uga_oda_5yr %>%
  mutate_if(is.numeric, round, 3)

crs_uga_oda_5yr <- crs_uga_oda_5yr %>%
  mutate_at(c("Avg. WASH Commitments Share", "Avg. WASH Disbursements Share",
    "Avg. BW Commitments Share of WASH", "Avg. BW Disbursements Share of WASH"),
    ~. * 100) %>%
  mutate_at(c("Avg. WASH Commitments Share", "Avg. WASH Disbursements Share",
    "Avg. BW Commitments Share of WASH", "Avg. BW Disbursements Share of WASH"),
    ~paste0(., "%"))

crs_uga_oda_5yr <- crs_uga_oda_5yr %>%
  mutate_all(~str_replace(., "NaN%", "-")) %>%
  mutate_all(~str_replace(., "NaN", "-"))

## Turn from wide to long

crs_uga_oda_5yr_l <- crs_uga_oda_5yr %>%
  pivot_longer(cols = starts_with("Avg"),
    names_to = "Variable",
    values_to = "Value")

crs_uga_oda_5yr_l <- crs_uga_oda_5yr_l %>%
  pivot_wider(names_from = "Period",
    values_from = "Value",
    id_cols = "Variable")

crs_uga_oda_5yr_l <- crs_uga_oda_5yr_l %>%
  rename("ODA Category" = "Variable")

## Export

stargazer(crs_uga_oda_5yr_l, summary = F,
          font.size = "small",
          column.sep.width = "2pt",
          rownames = F,
          title = "Total and WASH ODA Averages by Period (millions of 2019 USD)",
          label = "crs_uga_oda_5yr",
          type = "html",
          out = paste0(tables_directory, "crs_uga_oda_5yr.rtf"))

stargazer(crs_uga_oda_5yr_l, summary = F,
          font.size = "small",
          rownames = F,
          column.sep.width = "2pt",
          title = "Total and WASH ODA Averages by Period (millions of 2019 USD)",
          label = "crs_uga_oda_5yr",
          out = paste0(tables_directory, "crs_uga_oda_5yr.tex"))



